%--- main parameters
% SI System
%clear all;
function [p q] = Benchmark(option)
global U  h  g  p  q r  dt nt dx h1 x p_global tmax L steps y Q;
y0=[p;q;r];

dt=1;
I = eye(2*x,2*x);
Z1 = zeros(2*x,x);
Z2 = zeros(x,2*x);
Z3 = zeros(x,x);

M = [I Z1; Z2 Z3];
tspan = 0:dt:tmax;  % time vector, 
options1=odeset('RelTol',1e-6,'Stats','on'); 
options2=odeset('RelTol',1e-8,'Stats','on', 'Mass', M,'MassSingular', 'yes'); 

tic

if option ==1
    disp('ODE 45 Explicit Method Selected');
[t,y]=ode45(@waves,tspan,y0,options1);
else
    disp('ODE 15s Implicit Method Selected');
    [t,y]=ode15s(@waves,tspan,y0,options2);
    %[t,y]=ode15s(@upwinding,tspan,y0,options2);

end
toc


